#include "ErrorCheck.h"
#include <iomanip>
#include <cmath>



void printVec4(float *data, unsigned nData, std::ostream &out)
{
	for (unsigned i = 0; i < nData; i += 4)
		out << std::scientific << std::right << std::setprecision(2)
			<< std::setw(12) << data[i] << " "
			<< std::setw(12) << data[i+1] << " "
			<< std::setw(12) << data[i+2] << " "
			<< std::setw(12) << data[i+3] << std::endl;
}
void plasterBodies(float *data, unsigned nData, 
				float region, float mass)
{
	srand(0);
	for (unsigned i = 0; i < nData; i+=4) {
		//data[i] = (float) rand() / RAND_MAX;
		data[i] = ((float)rand() - (float)RAND_MAX / 2) / RAND_MAX * region;
		data[i+1] = ((float)rand() - (float)RAND_MAX / 2) / RAND_MAX * region;
		data[i+2] = ((float)rand() - (float)RAND_MAX / 2) / RAND_MAX * region;
		data[i+3] = mass;
	}
}
void describeError(float *aa, float *ab, unsigned nData)
{
	float dMax, dAvg;
	estimateError(aa, ab, nData, dMax, dAvg);
	std::cout
		<< "Error: max=" << dMax 
		<< ", avg=" << dAvg
		<< std::endl;
}
void estimateError(float *aa, float *ab, unsigned nData, float &dMax, float &dAvg)
{
	float d;
	double dSum = 0;
	unsigned iMax = 0;

	dMax = 0;
	for (unsigned i = 0; i < nData; i+=4) {
		if ((d = fabs(aa[i] - ab[i])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
		if ((d = fabs(aa[i+1] - ab[i+1])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
		if ((d = fabs(aa[i+2] - ab[i+2])) > dMax) {
			dMax = d;
			iMax = i;
		}
		dSum += d;
	}
	dAvg = (float)(dSum / nData);
}
void computeOnCPU(float *bodies, float *accels, 
				  float eps, unsigned nData)
{
	double r2, r3, rij[3];
	double tmp[4];


	for (unsigned i = 0; i < nData; i+=4) {
		//forces[i] = forces[i+1] = forces[i+2] = forces[i+3] = 0;
		tmp[0] = tmp[1] = tmp[2] = tmp[3] = 0.0f;
		for (unsigned j = 0; j < nData; j+=4) {
			rij[0] = bodies[j] - bodies[i];
			rij[1] = bodies[j+1] - bodies[i+1];
			rij[2] = bodies[j+2] - bodies[i+2];
			r2 =  rij[0] * rij[0]
				+ rij[1] * rij[1]
				+ rij[2] * rij[2]
				+ eps;
			r3 = r2 * sqrt(r2);
			tmp[0] += bodies[j+3] * rij[0] / r3;
			tmp[1] += bodies[j+3] * rij[1] / r3;
			tmp[2] += bodies[j+3] * rij[2] / r3;
		}
		accels[i] = (float)tmp[0];
		accels[i+1] = (float)tmp[1];
		accels[i+2] = (float)tmp[2];
		accels[i+3] = (float)tmp[3];
	}
}
